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Abstract. We derive a novel lattice Hamiltonian, the Molecular Hubbard Hamiltonian (MHH), which 
describes the essential many body physics of closed-shell ultracold heteronuclear molecules in their 
absolute ground state in a quasi-one-dimensional optical lattice. The MHH is explicitly time-dependent, 
making a dynamic generalization of the concept of quantum phase transitions necessary. Using the Time- 
Evolving Block Decimation (TEBD) algorithm to study entangled dynamics, we demonstrate that, in 
' the case of hard core bosonic molecules at half filling, the MHH exhibits an emergent time scale over 

which spatial entanglement grows, crystalline order appears, and oscillations between rotational states 
self-damp into an asymptotic superposition. We show that this time scale is a non-monotonic function 
of the physical parameters describing the lattice. We also point out that experimental mapping of the 
g . static phase boundaries of the MHH can be used to measure the molecular polarizability tensor. 

q \ 1. Introduction 



> 
00 



In recent years, ultracold atomic gases have provided near perfect realizations of condensed matter 
Hamiltonians, acting as quantum simulators [U [2] that allow the study of complex condensed matter 



^ . phenomena in a clean and highly controllable environment. Ultracold polar molecular gases, which have 

^ recently been brought to the edge of quantum degeneracy in their absolute ground state [HI Hj, offer 

^vq additional features over atomic gases, such as a large internal Hilbert space and a greater susceptibility 

^ to external fields via a permanent electric dipole. There have been a number of proposals on how to 

© use ultracold molecular gases for mimicking well-known Hamiltonians such as spin-1 lattice models [5] 
* * 

> Ultracold molecules have also been suggested as a model system for the study of strongly correlated 2D 
^ quantum phases [6] or for quantum information processing schemes [3, El [9]. However, these proposals 
frequently involve complex and yet-to-be implemented experimental techniques. In this article, we 
instead focus on the completely new quantum many body physics which results naturally from the 
simplest quantum lattice experiments that can be performed in the immediate future with established 
techniques in ultracold molecular quantum gases. 

Towards this end we derive a novel lattice Hamiltonian, which we refer to as the Molecular Hubbard 
Hamiltonian (MHH). The MHH describes the physics of an ultracold polar molecular gas in a ID optical 
lattice that is oriented using a DC electric field, giving rise to a resonant dipole-dipole interaction, and is 
driven between rotational levels using a microwave AC field. In particular, new aspects of our derivation 
include explicit dependence of hopping energy on the molecular polarizability tensor. This in turn allows 
a determination of the tensor elements, an outstanding experimental issue, from the borders of the static 
phase diagram of the MHH, which are identical to those of the extended Bose-Hubbard Hamiltonian [10] 
when a single molecular rotational level is occupied. 
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Beyond the statics, the MHH naturally has a dynamical component due to the AC driving fields, 
as well as an internal structure in terms of rotational modes which is inherently different from spinor 
atomic systems [TTI [T2] . We study this dynamical aspect with Time- Evolving Block Decimation 
(TEBD) p~3l [Hj, a newly developed entangled quantum dynamics algorithm which takes spatial 
entanglement (specifically, Schmidt number [15]) as a cut-off. We find an emergent time scale in the 
case of half-filling for hard core bosonic molecules. We emphasize that a quantum lattice model requires 
low filling (average number of particles per site), in contrast to a mean field lattice model, for which the 
filling would typically be quite high. Thus, although experiments can most easily access the mean field 
regime of hundreds of molecules per site with a single pair of counter-propagating laser beams, we look 
slightly ahead to the quantum regime, which will require two pairs of such beams in order to create an 
array of quasi- ID "tubes." A third pair is then used to create the lattice in each tube. This technique 
is already well established for ultracold atoms [16]. 

Dynamical aspects of quantum phase transitions are just beginning to be considered pEH HBj, and 
have so far been a limited area of study restricted to mean field considerations, due to lack of numerical 
tools. With the recent advent of entangled quantum dynamics algorithms, namely TEBD, dynamical 
properties of many-body systems are becoming amenable to numerical study. For example, TEBD has 
been used to address key questions such as the dynamics of a quantum quench [T9], [20] or the speed 
at which correlations propagate in a lattice [21]; these are not issues which can be studied with other 
dynamical methods such as dynamical mean field theory (DMFT) [22] . We give a brief review of TEBD 
in Sec. 3. The reader interested in computational details can find them in Ref. [23] . 

The first main contribution of this paper is to present a careful derivation of the Molecular Hubbard 
Hamiltonian. This is done in Sec. El with some previously known aspects of molecular physics relegated 
to | Appendix A , The second main contribution is to present an emergent time scale for half filling; 
although we treat the case of hard core bosons, the MHH can also be applied to fermionic molecules. 
To this end, in Sec. [3] we first give a brief explanation of TEBD and the quantum measures we use. 
Then, in Sec. [4] we present and analyze our simulations, with an accompanying convergence study 
in 



Appendix B[ Finally, in Sec. [5] we summarize. 



2. The Molecular Hubbard Hamiltonian 

The Molecular Hubbard Hamiltonian (MHH) is 

H = — ^2 tjJ'M Y {^i^J'M^iJM + h.C.) 
JJ'M 

+ Yl E JM ™iJM ~ * sin fat) ^ ftjM Yl (4/,M^J+l,M + h.C.) 
JM i JM i 

+ 2 Y U dd Y ^\jiM^iJ[M^\' j 2 M'<*i> J' 2 M' • (1) 

Jl, J[, J2, J'<2 (h^) 
M, M' 

where clum destroys a bosonic or fermionic molecule in the \£\ JM) state (defined below) on the i th 
lattice site, and the bracket notation (...) denotes that the sum is taken over nearest neighbors. The 
first term in Eq. (pD) corresponds to hopping both between sites and molecular rotational states with 
quantum numbers J, M. The second term represents the rotational energy along with rotational state- 
dependent energy differences due to a DC electric field. The third term corresponds to an AC electric 
field, making this a driven system. The fourth term corresponds to electric dipole-dipole interactions. 
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In the following subsections and |Appendix A| we justify Eq. (pQ) with a careful derivation and present 



the energy scales of each term. 



2.1. Derivation of the Molecular Hubbard Hamiltonian 



The full molecular Hamiltonian in second quantization is 




(2) 



The terms on the first line correspond to single-molecule effects: kinetic energy, rotation, the DC electric 
field which orients the dipole, the AC microwave field which drives transitions between rotational levels, 
and the far off-resonant optical lattice potential, respectively. The second line is the two-molecule 
resonant dipolar energy. The field operators ijj can be either bosonic or fermionic. We focus on the 
bosonic case for brevity. There are five key assumptions underlying our derivation, as follows. We 
consider all five assumptions to be reasonable for present and near- future experiments. 

(i) We consider ultracold closed-shell polar heteronuclear diatomic molecules, characterized by 
permanent dipole moment d and rotational constant B. The most experimentally relevant bosonic 
species in this category are SrO, RbCs, and LiCs [6]. The individual molecules are assumed to 
be in their electronic and vibrational ground states, and it is assumed that none of these degrees 
of freedom can be excited at the large intermolecular separations and low temperatures/relative 
energies that we consider. 

(ii) The molecule is assumed to have a 1 S ground state. The characteristic trapping potential length is 
chosen large enough compared to the internuclear axis to assume spherical symmetry, i.e. a locally 
constant potential. 

(iii) We neglect any intramolecular interactions (e.g., hyperfine structure), as they are typically very 
small for 1 S molecules [24]. 

(iv) We consider only the lowest three rotational levels. All AC fields will be sufficiently weak to allow 
this assumption. 

(v) We work in the "hard-core" limit where at most one molecule is allowed per site. This is enforced 
by strong dipole-dipole interactions on-site. We consider the lattice spacing large enough to include 
only nearest-neighbor dipole-dipole interactions. Other short-range interactions such as exchange 
or chemical reactions or long range interactions such as dispersion and quadrupole-quadrupole 
interactions are not considered. 

We proceed to follow the usual procedure [25] of expanding the field operators of our second- 
quantized Hamiltonian in a Wannier basis of single-molecule states centered at a particular discrete 
position iv 



where i is a site index and the sum is over all lattice sites. For our Wannier Basis we choose the single- 
molecule basis that diagonalizes the rotational and DC electric field Hamiltonians, spanned by kets 
\S\ JM). In this basis, which we refer to as the "dressed basis" (the DC field "dresses" the rotational 
basis) we have the field operator expansion 



(3) 



^2a iJM wjM (r-r<) = ^ d iJM \E\ JM)< . 



(4) 
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We note that such a basis, while highly efficient for the hard core limit we consider, becomes progressively 
worse for higher filling factors, till in the mean field limit the single-molecule basis, whether dressed or 
not, is so poor that many bands must be considered. Here we do not include a band index for simplicity, 
although the generalization of Eq. (pD) to include multiple bands is straightforward. 

This choice of Wannier basis associates the terms in Eq. (pD) to the terms in Eq. ([2j) as follows: 



* j,j',m = - J dr w k JM (r - i-) [H k[n + H opt ] Wj, M (r - r i+1 ) , (5) 

E JM = Jdr Wj M (r - r 4 ) [H rot + H DC ] ™jm (r - r^) , (6) 

-ttQjm sin (cut) = Jdr w k JM (r - i-) [H AC ] w J+1M (r ~ r i) , ( 7 ) 

U dd M,M ' = I drdr'wj, iM (r - r<) Wj^ M , (r ; - r m ) H dd (r - r') w JlM (r - r<) wj 2 m> (r 7 - r i+ i) ,(8) 

where the operators H\^ n) H opt , etc., are taken to be in position space representation. For the derivation 
of the single-molecule terms (rotational, DC electric field, and AC electric field) and discussion of the 
properties of our Wannier basis we refer the reader to | Appendix A , In the following sections we 
present the derivation of the tunneling (hopping) and dipole-dipole terms, which have new aspects not 
heretofore appearing in the literature [26] . 



2.2. Tunneling 

The tunneling term represents the sum of the molecular kinetic energy with the potential energy of 
the lattice. After expanding in the Wannier basis of Eq. (|4j), we find the effective tunneling Hamiltonian 

Hf = - t JJ'M Yl {^{j'M^i'.JM + h.c.) (9) 

where tj,j/,M was defined in Eq. (G3). To understand why this operator mixes states of different J, we note 
that the kinetic energy and (far off-resonant) optical lattice potential do not mix rotational eigenstates. 
Because our Wannier basis states are dressed and therefore superpositions of rotational eigenstates with 
different J, the tunneling operator in the dressed basis will mix J. Although the dressed basis makes 
the tunneling more complex to analyze, it simplifies other terms in the MHH, such as the DC term, and 
is in any case a more standard basis for analysis of the diatomic molecules we study here. Comparable 
basis changes are sometimes made in other quantum many body systems, where, for instance, particles 
and holes are mixed, or particles are paired. Note that, because we assume z-polarized fields, M is still 
a good quantum number. To discuss the actual form of the tunneling energies {tj,j/,M} we must first 
examine the interaction of a diatomic molecule with the optical lattice. 



2.3. Interaction with an Optical Lattice 

The charge redistribution that occurs when a molecule is subjected to a static, spatially uniform 
electric field E is reflected in its dipole moment d via the polarizability series 

dj = df^ + aj k E k + —PjkiEkEi + —Tj k i m E k EiE m + . . . (10) 

where the first, second, and third order coefficients aj k) (3j k u and Tj k i m are elements of the polarizability, 
hyperpolarizability, and second hyperpolarizability tensors, respectively. The polarizability tensor is a 
symmetric rank-two tensor with no more than six independent elements (less if molecular symmetry is 
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greater), and characterizes the lowest order dipole moment induced by an applied electric field. From 
this tensor we can form the scalar invariants 

-Trfi, (11) 



a 



3 
1 



(Aa)' = 2 [3Tr(« 2 ) 



(Tvaf 



(12) 



referred to as the polarizability and the polarizability anisotropy, respectively. Note that we use the 
tilde to clarify that a with elements ajk is a tensor, not a scalar - we reserve the accent circumflex (the 
"hat" symbol) for quantum operators. In linear molecules, such as diatomic molecules, the presence 
of only two distinct moments of inertia allows for the classification of a according to its components 
along and perpendicular to the internuclear axis, denoted a\\ and a±, respectively. In the presence of 
AC electric fields with frequency uj we speak of the dynamic polarizability tensor a (a;), with the series 
of Eq. (fTUl) being the zero frequency limit. The tensor a {uj) is, in general, complex, with the real part 
inducing a dipole moment and the imaginary part accounting for power absorption by the dipole and 
out-of-phase dipole oscillation. In the case of E diatomic molecules in their electronic and vibrational 
ground states [27] 



a{u) = a\\ {uj) e' ® e' + a± {uj) ^ (~l) A e 

A=±l 

where the e' q are molecule-fixed spherical basis vectors, 
polarizabilities are 

2 

d u X(v)-XZ(0) 



-A j 



(13) 



CL\\ 



= EE 



± u,v Ev?>{v) 



2 

4n(v)-xs(o) 



The parallel and perpendicular dynamic 



(14) 



zb v,v 



(15) 



respectively. In these expressions 4a(^)-xe(o) is the transition dipole moment from the ground state to 
the WY (v) state (following the usual diatomic molecular notation, A 6 {E, 11} = {0, 1} is the quantum 
number associated with the projection of the total electronic orbital angular momentum along the 
internuclear axis, i.e., in the molecule- fixed basis) and the sum over =p accounts for the near-resonant 
and typically far off-resonant terms. 

Transforming a from the molecule-fixed basis to the space-fixed basis using the transformation 
discussed in |Appendix A[ we find 



*>l) =£ E E ( 2 ^' + 1 ) 

P1P2 j=0,2 m=-j 

x Li (j + 2) (j - 1) - 4a ± 




1 



i)!(3 + i)! 



(16) 



where C$ is an unnormalized spherical harmonic, (. . .) denotes the Wigner 3-j coefficient [28J, and the 
are space- fixed spherical basis vectors. 
The interaction of the lattice with the molecule is represented by the Hamiltonian 

# opt (x) = -E: pt (r)-a'K).E opt (r) . 

If the electric field has polarization p in the space-fixed spherical basis then we find 



(17) 



#opt (x) 



lEopt (r)f 



(a ll +2a ± )C ( o }) + (-lf 



(l-p)!(l+p)! 



(a\\ — a±^j Cq 



(2) 



(18) 
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For light linearly polarized in the x-direction we obtain 

\2 



opt 



E pt (r) 
6 



2 (a\\ + 2a ± ) Cj 0) + (a,, - a±) (VsC. 



(2) 
2 



(19) 



whereas for light linearly polarized in the ^-direction we find 



^opt — 



|E opt (r) 
6 



-2 (a„ + 2a ± ) + (a„ - a ± ) (V6C<$ + 2Cf + V^Cf)] . (20) 



Since Cq^ = 1, these terms give a state-independent energy shift. The terms produce a tensor shift. 
Because the depth (in energy) of a typical optical lattice is much smaller than the energy of transitions 
between rotational levels (of order £>, as defined in Appendix A), we can ignore far off-resonant Raman 
coupling between different J manifolds and use only the diagonal matrix elements. The C\ term 

(2) 

and the C_ 2 will both mix M in the J > 2 manifolds, but do not affect the lowest two rotational 
levels, again, because we neglect Raman couplings. Thus x, y ) and z polarizations all have the same 

(2) 

Hamiltonian in this approximation. We can calculate the matrix elements of Cq in the field free basis 
using the Wigner-Eckart theorem to find 

|2 



(J'M'\H opt (r)\JM) = 



|Eopt W 



, a\\ 



+ (-if 



2ol) 

(a\\ — a± 



J(J+1) -3M 2 



(2J- 1) (2J + 3) 



Sjj'Smm'- (21) 



(l-p)\(l+p)\ 

In our effective Hamiltonian we choose right circular polarization for the z lattice, x polarization for the 
x lattice, and y polarization for the y lattice, where each "lattice" refers to a pair of counter-propagating 
laser beams used to create a standing wave. 

We consider the fields making up the optical lattice to have sinusoidal spatial profiles, resulting 
in sine-squared intensity profiles. In addition, we assume that the y and z lattices are tight, meaning 
that the molecules are strongly confined at the potential minimum (for a red-detuned trap). This tight 
confinement allows us to approximate them via a Taylor series, e.g., sin 2 (k z z) ~ k 2 z z 2 in the vicinity of 
the molecule. Using the above results, the matrix elements of the Hamiltonian for the optical lattice 
can be written 



(J'M'\H opt (v)\JM) 



|E op t (y)| 2 k 2 y y 2 + |E opt (x)| 2 sin 2 (k x x) 



a + 2Aa 



J(J+ 1) -3M 2 
(2J- 1) (2J + 3) 



|E op t (z)| 2 ^ 2 



a 



Aa 



J(J+1)-3M 2 
(2J- l)(2J + 3) 



S.TJtS 



JJ'VMM' 



<5 J J> 'MM' 



(22) 



or, more compactly, as 

(J'M'\H ovt (r) \JM) = \-a% |E opt (y)| 2 k 2 y y 2 - a% |E opt (x)| 2 sin 2 (k x x) 



Sjj'S 



J J' VMM' 



- |E op t (z)| 2 ay M k 2 z z 2 5j Jf 5 M M> 



by defining 



a 



W = 

JM — 



a + 2Aa 



a 



W = 

JM — 



a 



Aa 



J(J+1)-3M 2 
(2J- l)(2J + 3) 
J(J+1)-3M 2 ' 
(2J- l)(2J + 3) 



We now define, as is customary, the "lattice heights" in the x, y, and z directions, respectively, as 



= - |E opt (x 



a 



(*) 

JM") 



(23) 

(24) 
(25) 

(26) 
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v y w= -|E opt (y)| 2 41, 



y(.JM) _ 



IE. 



(27) 

. opt (z)r4t. (28) 

The tight confinement in the transverse (y and z) directions strongly suppresses tunneling in these 
directions, making the overall lattice effectively ID along x. 

From Eqs. (I24l) - (l25l) . it is apparent that different rotational levels experience different trapping 
frequencies and different tunneling energies. To make this clearer, we parse our full field-free tunneling 
matrix element as 



t 



JM 



I 



drw* JM (r 



dr w 



JM 



r — r 



+ J drw 



JM 



r») [H kin + H opt ] w JM (r - r i+1 ) 
-H kin + K (JM) sin 2 (k x x 2 )] w JM (r - r i+1 ) 

r - r< ) [vj jM ^y + vy^ky] WJM ( r - r<+1 ) 



Defining 



tf M = J dr w* JM (r - n) [-//kin + Vj JM) sin 2 (k x x 2 )] w JM (r - r i+1 ) , 



t 



(trans) 

JM — 



drw* JM (T-Ti) 



2^2 



vy^k-z 



2„2 



w JM (r - r i+1 ) , 



(29) 
(30) 



we proceed to compute each piece separately. 

In the evaluation of the first integral, Eq. (|29|) we assume that the Bloch function of a molecule in 
the sinusoidal optical lattice is a Mathieu function along x. This may seem to contradict our assumption 
of spherical symmetry in the above derivation. However, the assumption of spherical symmetry (i.e. 
a locally constant potential) need only hold on the order of an internuclear axis (~ 5 A) near the 



molecule. In contrast, on the order of the characteristic lattice length ^Jh/fiuj opt the rigid-rotor molecule 
is indistinguishable from a point particle (such as an alkali atom), and so spherical symmetry is not 
required. With this understanding, we recognize tjj^ as the expression for the hopping energy for point 
particles in optical lattices [29] with the additional feature that the lattice height along the quasi- ID 
direction Vq = Vy M ^ is dependent on J through the polarizability tensor. Thus, altering the expression 
from the theory of point particles in optical lattices, we obtain the result 

/ | jr(JM) 



t 



(0) 
JM 



/vy M >\ 



B 



exp 



-C^ 



V 



Er \ E R 

where A = 1.397, B = 1.051, C = 2.121, and 

E R = h 2 k 2 x /2m 
is the recoil energy. 



\ Er 



(31) 



(32) 



For the second integral, Eq. (I30l) . we approximate the Wannier functions with the ground state of 
a simple harmonic oscillator 



w(y) 



-1/2 



7T 1 ^ 4 exp 



-y 



(C) 
(C) 

where the harmonic oscillator lengths are given by 



-1/2 



W [ZJ 



7T 1 ^ 4 exp 



,2 (e> 



(33) 
(34) 
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1 T 71 /T\ 

\JM) 




o W / - 


1 r\r\\ 

|00) 


1 


1 


1 10) 


1.715 


0.642 


1 -i i -i \ 

|1±1) 


0.642 


1.178 


1 on\ 

|20) 


1.511 


f\ <~7 A A 

0.744 


|2± 1) 


1.255 


0.872 


|2±2> 


0.488 


1.255 



Table 1. Values of the polarizabilities for LiCs in different rotational states | JM). 
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\E opt \ 2 a/3E R 



Figure 1. Dependence of the field- free tunneling (hopping) coefficient on rotational state and lattice 
height. 



Then 



, (trans) 

t jM ex P 



A 2 



+ 



a 



a 



JM 

(*) 
JM 



exp 



A 2 



(36) 



where A is the wavelength of the optical lattice. Because we consider tight traps such that the lattice 
height in the y and z directions is much greater than the lattice height in the x direction, V y ~ V z ^> V Xl 
this contribution is exponentially suppressed compared to t^j M , and so we neglect it. Thus, 

B ( 



t 



JM 



t 



(0) 
JM 



E 



R 



E, 



A 



IE, 



opt I 



a 



JM 



R 



E 



R 



exp 



-C 



\ 



\ 



IE, 



opt | 



a 



JM 



E 



R 



(37) 



This is equivalent to the array of tubes we discussed in Sec. [H where each tube is isolated from its 
neighbors. 

Using tabulated values of the polarizabilities for LiCs[30j as given in Table [U we find that, for a 
reasonable lattice height V(°°'/Er ~ 10, the tunneling term for the |11) state is only about 20% of that 
in the 1 00) state, as shown in Fig.[U For LiCs in a red-detuned optical lattice of wavelength A = 985nm, 
E R = 2tt x 1.46ft kHZ. Typical values of the lattice heights are V x - IQE R) V y) V z - 25E R [33]. 

We reiterate that the above matrix elements and tunneling energies {^jm} have been computed in 
the field-free basis for simplicity. To transform to the dressed basis, we use the unitary matrix with 
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dressed eigenvectors as columns, recovering Eq. ([9]), where the tunneling matrix element is no longer 
diagonal in J. 

2.4- Dipole-Dipole Interactions 

The induced dipoles from the DC field give rise to a resonant dipole-dipole interaction. The 
Hamiltonian for this interaction in the two-site dressed basis spanned by \S ; J\M\J 2 M 2 ) is 

Hdd = ~ U dd Y^iJxM^iJ'^M^i'j^M'^i'J^M', (38) 

Ji, J[, j 2 ,j' 2 (i,i>) 

M, M' 

where we have defined 

U dd M,M ' = J drdr'wj^ M (r - r^) Wj^ M , (r' - r i + 1 ) H dd (r - r') w JlM (r - r^) wj 2 m> (r' - r i+ i) , (39) 

and for notational simplicity we have suppressed the S subscripts. Note that because of our choice of 
polarizations of the optical lattice and AC and DC electric fields, Mi = M 2 = M and M[ = M f 2 = M f . 

The resonant dipole-dipole interaction between two permanent dipoles d\ and d 2 whose respective 
centers of mass are separated by a vector R in the space-fixed frame is 

di • d 2 - 3 (e R ■ dj (d 2 ■ e R ) 

H dd = ^ , (4U) 

where e# is a unit vector in the direction of R. Using standard angular momentum recoupling we recast 
this in spherical tensor notation as 



(2) 



H dd = -TxfE ( R ) di ® d 2 j , (41) 

where (T)j^ denotes the component of the rank-fc spherical tensor T that has projection q along R, 
C$ (R) is an unnormalized spherical harmonic in the polar coordinates defined with respect to R, and 
we have defined the tensor product of the vector operators di and d 2 as 

d 1 ®d 2 l (fc) = ^(l,m,l, (? -m|A: (? )(d 1 ) (1) (d 2 ) (1) . (42) 

m 

In the last line, (j'i,toi, j 2) m 2 \J ) M) is a Clebsch-Gordan coefficient. We now take matrix elements of 
Eq. ( 14B in the two dressed-molecule basis \S] J\M\, J 2 M 2 ) ) where molecule 1 is on site i and molecule 
2 is on site i + 1, yielding 

(£; J[M[, J' 2 M' 2 \H dd \£- J,M U J 2 M 2 ) = £ (-1)" C<$ (R) 

\E-J l M l ){£-J' 2 M' 2 \(d 2 f ) 

Because our DC field is polarized along z, only (di)o^ and (d2)o matrix elements are nonzero, enforcing 
/i = 0, m — 0. With this in mind, the interaction takes the particularly simple form 



x 



£<1, m, - m|2/x)<£; J(M(| (d^ \S; JiM^f; J f 2 M' 2 \ (d 2 ) (1) |£; J 2 M 2 ) . (43) 



(£; J[M[, J' 2 M' 2 \H dd \£- J 1 M 1 , J 2 M 2 ) = -]|cf (R) 

x <1,0,1,0|20)<£; J(M!| (dx)^ |£; JiMx}(£; J^M 2 | (d 2 )^ |£; J 2 M 2 ) (44) 
= (£; J(Mx| (dx)^ If; JiM^S; J^M 2 | (d 2 )^ |£; J 2 M 2 ) (— • (45) 
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The intermolecular axis plays a crucial role in the sign of the interaction. Two molecules oriented along 
the intermolecular axis attract if their dipoles are parallel and repel if their dipoles are antiparallel. Two 
molecules oriented perpendicular to the intermolecular axis, on the other hand, repel if their dipoles are 
parallel and attract if their dipoles are antiparallel. The DC field that orients the molecules in our setup 
is polarized along z, perpendicular to the intermolecular quasi-lD axis x. This gives rise to repulsive 
interactions for positive dipole matrix elements. With this geometry the dipole potential becomes 

(£; J[M[, J' 2 M' 2 \H dd \£- J,M ± , J 2 M 2 ) = 

^3<5; J[M,\ (dx)^ \S; J^M^E ; J' 2 M 2 \ (d 2 )f \S; J 2 M 2 ) , (46) 

yielding 

v = 8_ {£ _ J[m ^(D {£ _ JiMi){s . j, M2{ ^« ]e . hm j (47) 

where A is the wavelength of the optical lattice. 



2.5. Energy Scales 

We proceed to clarify the energy scales associated with each term in Eq. (GQ). Between previous 
discussion in Sec. [2] and that of | Appendix A[ all terms in Eq. (OQ) are now clearly defined. The energy 



scales of the dressed basis are £>, the rotational constant, which is roughly 60ft GHz, and dEvc, which 
is of order 1 — 10£>. The DC term has no length scale associated with it because the field is uniform, 
and the length scale of the rotational term is the internuclear separation, on the order of angstroms. 
The relative contribution of the DC electric field and rotational terms in Eq. (pp) are expressed through 
the dimensionless parameter 

Pdc = d£ DC /B, (48) 

the ratio of the DC field energy to the rotational level splitting. 

The energy scales of the AC term are Tiuj, where u is the angular frequency of the driving field, and 
dS^c- The scale huj is of order 2B for small /?dc ^ 1? and of order £>V/?dc for large /?dc 3> 1- The AC 
field energy d8^c is of order 0.5huu. The single-molecule time scale associated with rf^AC is the Rabi 
period, the time it takes for the population of a two-level system to cycle once, as seen in Figure 3(a)[ 



In real time, this is on the order of lOps for the parameters in the preceding paragraph. The time scale 
associated with cu is the time scale on which the small oscillations in Figure 3(a)| occur, of order 0.5ps. 



The length scale of the AC field is on the order of centimeters, and so we can neglect this in light of 
the micron length scale of the trap. 

The tunneling term has several scales. The optical lattice near the point of confinement has a 
length scale given by the harmonic oscillator length /(^ ~100nm and an energy scale of Er ~1ATi 
kHz. The energy scales of the tunneling operator proper are given by the {^jj/m} which are of order 
10~ 1 -10 _2 £'fl ~100ft Hz for the given recoil energy. 

There are also many scales for the dipole term. For the B and d specified in the first paragraph 
of this section and /?dc = 1.9, the characteristic length scale where the dipole-dipole energy becomes 
comparable to the rotational energy is 



r B 



2 \| 

(£;00|d|£;00) /B) , (49) 



approximately 348 Bohr radii (18.4nm). Outside this region the Born-Oppenheimer adiabatic 
approximation is easily fulfilled [6]. Since the length scale of our optical lattice is of order /xm, we 
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Term 


Length scale 


Energy scale 


Rotation 


internuclear distance ~ 1 A 


B ~ 60ft GHz « 2cm" 1 


DC field 


N/A, uniform 


d£ DC ~ 120ft GHz « 4cm" 1 


AC field 


2ttc/uj ~ 1cm 


ftw - 30ft GHz « 1cm" 1 


Kinetic 


CI ~ lOOnm 


£ R - 1.46ft kHz 


Tunneling 


Lattice spacing^ 1/im 


{^j'jm} ~ 100ft Hz 


Resonant Dipole-Dipole 


energy comparable to B 
at tb — 348 Bohr radii 


(£*; 00|d|£:; 00) / (Umi) 3 ~ 1.2ft kHz 
for nearest neighbors 



Table 2. Comparison of energy and length scales for the Molecular Hubbard Hamiltonian of Eq. (pQ). 



are justified in working within the Born-Oppenheimer framework. For the same parameters, the length 
scale where the off-resonant van der Waals potential Cq/t 6 « — d 4 /(6£>r 6 ) becomes comparable to the 
dipole-dipole interaction is 

r v dw = (2 \C 6 \ I \(£- 00|d|£; 00>| 2 )* . (50) 

This length is very small, on the order of tens to hundreds of Bohr radii. Outside of this region the 
resonant dipole potential dominates and the intermolecular force is repulsive. This repulsion enforces 

2 

the hard-core limit. The energy scale of the dipole-dipole force is (£; 00|d|£; 00) /A 3 ^ 1.2ft kHz, with 
higher J being an order of magnitude or so lower for small /?dc? an d of the same order for large /?dc 
(see Fig. 1(a)). 

To summarize, the scales of the problem are shown in Table El 



2.6. Novel Features of the Molecular Hubbard Hamiltonian 

The MHH, Eq. (pQ), has a number of novel features which distinguish it from the Hamiltonians 
typically considered in the quantum lattice and condensed matter literature (32], [25] . First, the tunneling 
energies {ij,j'M} not only depend on the rotational level J, M but even change rotational states from J 
to J f . This is due both to the polarizability tensor's dependence on rotational level, and to the dressed 
basis. This differs from other Hubbard models which consider spin degrees of freedom, as tunneling does 
not occur between spin states - hopping does not cause spin transitions. If we consider populating a 
single mode (e.g. J = 0, M = 0) in the Q — > limit, then Eq. (GD) becomes the extended Bose-Hubbard 
Hamiltonian, and the phase diagram is known [33lfT0] . This gives ideas of how to characterize the static 
phases of the MHH. However, because the tunneling energy depends on J, the borders of the phase 
diagram will depend on the rotational state of the system. We will discuss this property and provide 
an application in Sec. [H 

Second, the Hamiltonian is fundamentally time-dependent because it is a driven system. This 
allows for the study of dynamic quantum phases, requiring the concept of a quantum phase diagram to 
be generalized to an inherently time-dependent picture. In a case study for hard core bosonic molecules 
at half filling presented in Sec. [H we show that the MHH has an emergent time scale. 
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3. Methods 

3.1. Time-Evolving Block Decimation 

The Time-evolving Block Decimation algorithm (TEBD) is a new method [32 [35] designed to 
study the dynamics of entangled quantum systems. The essential idea of TEBD is to provide a moving 
"spotlight" in Hilbert space which tracks a dynamical system. The portion of the Hilbert space so 
illuminated is an exponentially small fraction of the full Hilbert space; this is justified by the fact that 
real, physical quantum many body systems, especially in real materials, typically explore only a small, 
lowly-entangled part of the total Hilbert space. 

In fact, TEBD moves the full quantum many-body problem from the NP-complete complexity class 
to the P class through an exponential reduction in the number of parameters needed to represent the 
many body state. We can understand the possibility of this reduction through an analogy to image 
compression. Present digital cameras are capable of producing a roughly 3000 x 3000 array of pixels. 
Downloading the images from such a camera, one notices that there are far less than 10 Megapixels 
worth of data per image. Image compression algorithms such as JPEG produce images of remarkable 
quality with only a small fraction of the raw data. The reason that these algorithms are so effective 
is that a physical image, as opposed to a random 2D pixel array, is not the "most common" or most 
probable image; it contains a great deal of structure and regularity. In the same way, physical states 
in Hilbert space tend to be lowly entangled (by some entanglement measure), even though a general 
state in Hilbert space has a much larger probability of being highly entangled. There is no general 
proof of this fact, just as there is no guarantee that an image will come out perfectly crisp after JPEG 
compression; it is simply a trend observed in many-body quantum systems. 

To be slightly more specific, TEBD performs a partial trace over a particular bipartite splitting of 
the lattice, and then keeps the x largest eigenvalues of the resulting reduced density matrix. The cut-off 
parameter \ 18 based on the Schmidt measure [15], and so it also serves as a measure of the degree of 
spatial entanglement. This idea is not unique to TEBD. In fact, the density matrix renormalization 
group (DMRG) method first proposed by White [36] did something analogous years before. TEBD's 
innovation is that at each time step it re-optimizes the truncated basis (thus the "moving spotlight"). 
The Schmidt number is just the number of non-zero eigenvalues in the reduced density matrix, and so is 
an entanglement measure natural to quantum many body systems. The parameter x is the number of 
non-zero eigenvalues in the reduced density matrix that TEBD retains. It is the principal convergence 
parameter of the algorithm, both in entanglement and in time. Although the time-propagation method 
we use is Trotter-Suzuki [37], it turns out that, due to a normalization drift, x controls convergence at 
long times. 

With x interpreted as an entanglement measure, we can say that TEBD treats the system not as 
a wavefunction in a d L -dimensional Hilbert space (L is the number of lattice sites), but as a collection 
of wavefunctions in d 2 -dimensional two-site spaces that are weakly entangled with the environment 
created by the rest of the system. To facilitate this viewpoint, we replace the d L coefficients of the 
full many-body wavefunction with L sets of (dx 2 + x) coefficients corresponding to the wavefunctions 
of each bipartite splitting. The most computationally expensive portion of the TEBD algorithm is 
typically the diagonalization of these local coefficient matrices at a cost of O (d 3 x 3 )- Looping over all 
L — 1 bipartite splittings and evolving the system for a total time tf in time steps of length St, one 



This scaling can be greatly improved by the presence of conserved quantities. When a conserved 
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quantity exists in the system we are able to diagonalize reduced density matrices corresponding to 
distinct values of this conserved quantity independently, which can result in significantly smaller reduced 
density matrices to diagonalize. Implementing this idea, scalings of O (x 2 ) have been reported for fixed 
d [38]. In addition, conserved quantities in the presence of selection rules can reduce the local dimension. 
For example, in the case of the MHH, z-polarized electric fields disallow transitions from a particular 
M to any other. If we begin with all molecules in a particular M state, this allows us to restrict our 
attention only to states with this M. In our numerics we conserve both the projection M, and the total 
number N. Furthermore, to match our hard core requirement, we allow only zero or one molecules per 
site, so that the local dimension is d < R+ 1, R being the magnitude of the greatest angular momentum 
that we consider (note that the local dimension rf, mentioned only here in Sec. 13.11 bears no relation to 
the permanent electric dipole moment d used throughout the rest of our treatment). 

A more detailed description of TEBD can be found in Ref. [23]. We also recommend Ref. [39], 
besides Vidal's original papers [341 [35]. 



3.2. Quantum Measures 

We use a suite of quantum measures to characterize the reduced MHH, Eq. f |55l ) below. The few- 
body measures we use are (nf), the number in the J th rotational state on the i th site, E = (H), the 
expectation of the energy, and j j {n J ) ) the average number in the J th rotational state per site (L is the 
number of lattice sites). The latter is a J-dependent filling factor. The many body measures we use 
include the density-density correlation between rotational modes J\ and J 2 evaluated at the middle site 

AJ1J2) ( \ L 1 A _ /-{Ji)AJ2)\ 



9 [ 2 1J2) [l^i) = (n[ J J2) ) ~ (n[p(nr), (51) 

where lq\ is the floor function, defined as the greatest integer less than or equal to q. As an entanglement 
measure we use the Meyer Q-measure (40], SH [42] 



Q = 



d 



d-l 



M 

l-£Tr(/3 (fc) )' 



(52) 

k = l 

where p^ is the single-site density matrix obtained by tracing over all but the k th lattice site, and the 
factor outside of the bracket is a normalization factor id is the on-site dimension) . This gives an average 
measure of the entanglement of a single site with the rest of the system. The Q-measure can also be 
interpreted as the average local impurity (recall that the Tr(p 2 ) = 1 if and only if p is a pure state). 

To determine what measures we can use to ascertain the static phases of our model we reason by 
analogy with the extended Bose-Hubbard Hamiltonian where we know that the possible static phases 
are charge density wave, superfluid, supersolid, and Bose metal [10]. The charge density wave is an 
insulating phase appearing at half integer fillings which has a wavelength of two sites. Like the Mott 
insulating phase, it has an excitation gap and is incompressible. While the extended Bose-Hubbard 
Hamiltonian has only one charge density wave phase due to the presence of only one species, the MHH 
has the possibility of admitting several charge density wave phases due to the presence of multiple 
rotational states. As such, we define the structure factor 

^ lJ2) = ^E(-l) |l " jl (^ (Jl) ^ J2) ), (53) 

where TV is the total number of molecules. We recognize this object as the spatial Fourier transform 
of the equal-time density- density correlation function between rotational states J\ and J 2 , evaluated at 
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the edge of the Brillouin zone. This measure is of experimental interest because it is proportional to 
the intensity in many scattering experiments, e.g. neutron scattering [13]. Crystalline order between 
rotational states J\ and J 2 is characterized by a nonzero structure factor S^ JlJ2 \ The charge density wave 
is the phase with crystalline order but no off-diagonal long-range order as quantified by the superfluid 
stiffness of rotational state J 

(note that p s bears no relation to the density matrix p) . If both the structure factor and the superfluid 
stiffness are nonzero, the phase is called supersolid. If both the structure factor and the superfluid 
stiffness are zero, the phase is called Bose Metal. Finally, if the structure factor is zero and the 
superfluid stiffness is nonzero, the phase is superfluid. In one dimension the entire superfluid phase is 
critical, and so there is no order parameter [10]. 



4. Case Study: Hard Core Bosonic Molecules at Half Filling 



In the following, we consider a particular case of Eq. (pD) for dynamical study. We choose the hard 
core case, which can occur naturally due to strong on-site dipole-dipole interactions, and half filling, 
which is an interesting point in a number of models, including the repulsive Fermi-Hubbard Hamiltonian 
and the extended Bose-Hubbard Hamiltonian discussed in Sec. 13.21 For example, in the latter case, the 
charge-density- wave phase requires a minimum of half-filling [10] . 

If we assume that our system begins in its ground state (J = 0, M = 0) we need only include 
states which have a dipole coupling to this state. For z-polarized DC and AC fields, this means we only 
consider M = states, yielding the reduced Hamiltonian 

ij + h.c.) + e j Y - n sin M) Y q j Y + h - c -) 



H = -J2tjj> Y (4,j' a * 



i 



(55) 



rr /i,J(,J 2 ,J 2 -t - -t - 
U dd 1^ a iJ 1 a U , 1 a i'J 2 a i , J^ ■ 

<M') 

This is the specific case of the MHH that we study using TEBD. 

A matter of practical concern, as apparent in Table El is the large disparity between the timescales 
of the first three (Rotational, DC, and AC) and the last three (kinetic, tunneling, and Dipole-Dipole) 
terms. The accumulation of error resulting from truncating the Hilbert space at each TEBD timestep 
causes the algorithm to eventually fail after a certain "runaway time," making studies over long times 
intractable [Sj. This invites a multiscale approach in the future [351 H6j. In our current numerics we 
artificially increase the recoil energy and dipole-dipole potential to be of the order of the rotational 
constant in order to study Eq. ( 1551 ) using TEBD. In particular, we take 



10B 



(£;J[\d\£; ■/!>(£; J£|d|£; J 2 > , 



(56) 



105 



x exp 



V 1 + 2 



Aa J(J+1) 
a (2J+l)(2J + 3) y 



1.051 



-2.121 



\ 



V 1 + 2 



Aa J(J + 1) 
a (2J+l)(2J + 3) 



(57) 
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where the dimensionless variable r] becomes an ersatz "lattice height." To see the scaling more explicitly, 
we compare the above with the actual expressions for the MHH parameters 



dd 



(£• J[\d\E; J x ) (E; J' 2 \d\E; J 2 ) 

(E-J^E-J^iS-J'^e-h)/^ 
Aa J(J+1) 



A 3 

( 2mE R d 4 / 3 



1 + 2- 




3M 2 



a (2J-l)(2J + 3) 



1.051 



x exp 



1 + 2 



Aa J(J+ 1) -3M 2 
a (2J-l)(2J + 3) 



(58) 
(59) 

(60) 
(61) 



2mE R d 4 / 3 / (ft 2 * 2 )] 5 = 105 for this E R , we 



If we now scale E R to be 105/1.397 and set d such that 
recover Eqs. (1551) and (1571) provided we make the definition 

r, ee - |E opt (x)| 2 a/ (3E R ) = V^a/ (3E R a%) . (62) 

Since this dimensionless parameter plays the same role as the quasi- ID lattice height scaled to the recoil 
energy did in the actual MHH, we refer to it as the lattice height. For the polarizability tensor, we 
choose Aa/a = 165.8/237, corresponding to LiCs [30]. This rescaling does not change the qualitative 
static and dynamical features of Eq. f |55l ): it only makes Eq. f |55l ) treatable directly by TEBD, without 
multiscale methods. 

First, we point out that if we consider populating a single rotational state (e.g. J = 0, M = 0) 
in the Q — > limit, then Eq. f |55l ) becomes the extended Bose-Hubbard Hamiltonian, and the phase 
diagram is known [33l [TO] . Because the tunneling energy is different for different rotational states (see 
Eq. (EIJ)) and this difference depends only on the properties of the polarizability tensor, we can relate the 
borders of the phase diagram for different rotational states to properties of the polarizability tensor. The 
MHH thus gives a means to measure the polarizability tensor, a standing issue in experiments [47]. Our 
calculations in Sec. [2] can be used to compare directly to the phase diagram from the literature [33l fTO] . 
In fact, this aspect of our work, unlike the simulations below, is not restricted to ID. 

However, our main focus at present is on the dynamics of the MHH. In the following numerical 
study, we explore dynamics as a function of the physical characteristics of the lattice, namely, number 
of sites L and effective lattice height rj. Specifically, we study L = 9, 10, and 21 lattice sites with 
iV=4, 5, and 10 molecules, respectively, and rj ranging from 1 to 10. We fix the dipole-dipole term as 
in Eq. f |56l ). and fix the DC field parameter to be /?dc — 1-9. While /?dc = 1-9 may not correspond to a 
physically realizable situation, its exploration provides insight into the MHH. 

The Rabi oscillations between the J = and the J = 1 states damp out exponentially in the 
rotational time t r = Bt/Ti as 

(h ) = a -b e~ tr/r cos (c t r ) , (63) 
(ni) = a 1 -b 1 e~ tr/T cos (dt r ) , (64) 

with some characteristic time scale r, as seen in Fig. El We note that an exponential fit has a lower 
reduced chi-squared than a power-law, or algebraic fit. We also tried fit functions where the oscillations 
do not decay to zero, but rather persist with some asymptotic nonzero amplitude. We find that the fit 
functions Eqs. f |63l ) and f |64l ) above fit the data better as quantified by the convergence properties of the 
algorithms used, as discussed in|Appendix B, 
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L 9,7/ 1,/?dc= 1.9 

— (hooj/L 

— (h 1Q )/L 




(a) Site-averaged population vs. rotational time for 9 sites. 
Note the general theme; a gradual decrease (increase) of the 
maxima (minima) of oscillations. 



L= 10,/; = l,/3 D c= 1.9 



— (hooj/L 

— (h 1Q )/L 




20 40 60 80 100 120 140 160 



Bt/Ti 

(c) Site-averaged population vs. rotational time for 10 sites. 
Note that there is no significant difference between an odd 
and even number of sites. 

L = 21,77 = l,/3 D c= 1.9 

— (hooj/L 

— (h 1 o)/L 




Bt/Ti 



(e) Site-averaged population vs. rotational time for 21 sites. 
Note that there is no significant difference between this and 
the smaller system sizes. 



16 




Tiw/B 

(b) Squared modulus of Fourier transform of 
site- averaged J = population vs. rotationally 
scaled frequency for L = 9 sites. The arrow 
denotes the Rabi frequency Ooo- 

L = 10,77 = i,Pbc= 1-9 
10 | 




Tiw/B 

(d) Squared modulus of Fourier transform of 
site- averaged J = population vs. rotationally 
scaled frequency for L = 10 sites. 

L = 21, rj = l,/3 D c= 1-9 

10 




10 -io| ^ ^ , 1 

0.5 1 1.5 2 

hjj/B 

(f) Squared modulus of Fourier transform of site- 
averaged J = population vs. rotationally scaled 
frequency for L = 21 sites. 



Figure 2. Dependence of site-averaged number on lattice size L. For this set of parameters, the site- 
averaged J = and J = 1 populations appear to asymptotically approach quarter filling. The J = 2 
mode is populated slightly by off resonant AC couplings. The peak near the left side of the Fourier 
transform plots is the Rabi frequency £7qo, denoted by an arrow. 
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L = 9,7] = l,/5 D c= 1.9 



9,?/= 1,/?dc= 1.9 




20 40 60 80 100 120 140 160 

Bt/% 

(a) Structure factors vs. rotational time for 9 sites. Note the 



similar asymptotic behavior to the populations in Fig. 2(a) 



L = 10,7]= i,/? DC = 1.9 




20 40 60 80 100 120 140 160 

Bt/% 

(c) Structure factors vs. rotational time for 10 sites. There is 



no significant difference in the and S^ L1) between even 

and odd L. For the difference in see Fig. |3(f)} 



L = 21,71= 1 5 /?dc= 1-9 




20 40 60 80 100 120 140 160 

Bt/% 

(e) Structure factors vs. rotational time for 21 sites. Note 
the lack of significant difference with the smaller odd system 




(b) Squared modulus of Fourier transform of 
vs. rotationally scaled frequency for L = 9 
sites. Note the similarity with Fig. |2(b)| above. 



9,7]= 1,/?dc= 1-9 




(d) Squared modulus of Fourier transform of 
^i 10 ^ vs. rotationally scaled frequency for L = 9 
sites. Note the absence of the Rabi frequency. 

V = hPr>c= 1.9 




100 

Bt/% 

(f) Comparison of the correlation structure 
factor for odd and even numbers of sites. Note 
that the even site (exactly half filling) structure 
factor grows faster and larger than the odd site 
(slightly less than half filling) structure factor. 



Figure 3. Dependence of structure factors within and between rotational states J on the number of 
lattice sites. We do not consider the off-resonant J = 2 and higher rotational states because they have a 
very small occupation; J = 2 is shown explicitly in Fig. O 
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L = 21,r/ = 5,/? DC =1.9 




21,r/ = 5,/? DC = 1.9 



Bt/Ti 



(a) Site- aver aged population vs. rotational time for 21 
sites with r] = 5. Note that the J = and J = 1 states 
now appear to converge to different fillings. 



I = 21,7] = io,/? DC = 1.9 




Bt/Ti 



(c) Site-averaged population vs. rotational time for 
21 sites with n = 10. Note the similarity to the 
n = 1 case (Fig. 2(e)[ ) and the difference from the 
rj = 5 case(Fig. 4(a) )-the asymptotic behavior is not a 
monotonic function of the lattice height. 




(b) Squared modulus of Fourier transform of 
(noo) vs. rotationally scaled frequency for L = 21 
sites and n = 5. Note the presence of several 
new frequencies not observed in the n = 1 case 
(Fig. |2(f)| ). In particular, Qqo, 2^oo, and 3^oo, 
are denoted by arrows. 



21,7? = 10, /3 DC = 




(d) Squared modulus of Fourier transform of 
(noo) vs - rotationally scaled frequency for L = 21 
sites and n = 10. Note that the frequencies that 
emerged during n = 5 have persisted. 



Figure 4. Dependence of the asymptotic behavior of rotational state populations on the lattice height 
n. 



The time scale r also describes the decay of physically measurable quantities, for example the 
structure factors as defined in Eq. (|53l) and illustrated in Fig. We show the emergent time scale r 
for various lattice heights and systems sizes in Table [31 

Examining Fig. El one observes that the driven system approaches a dynamical equilibrium that 
is a mixture of rotational levels. The time scale with which the system relaxes to this equilibrium, 
r, cannot be determined from the single-molecule physics, and so we refer to r as an emergent time 
scale. For the low lattice height rj — 1, the populations of the first two rotational states appear to 
oscillate around and asymptotically converge to roughly quarter filling, with J = 1 being lower due 
to contributing to population of J = 2 via an off-resonant AC coupling (Fig. 2(a)). For r] = 5, the 
asymptotic equilibrium is an uneven mixture of rotational states that favors occupation of the J = 
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L = 21,j7 = 5,/?dc=1.' 




L = 21,/3 DC = 1.9 



Bt/h 



(a) Structure factors vs. rotational time for 21 sites 
with T] = 5. 

L = 21, 77 = 10,/3 D c = 1.9 




Bt/h 

(c) Structure factors vs. rotational time for 21 sites 
with 77 = 10. Note the similarity of 5i° ^ and aS^ 11 ^ to 



the 77 = 1 case (Fig. |3(e) ). Note also that is now 
nonzero, and is periodic with the Rabi frequency ^00 
at short times and twice the Rabi frequency at long 
times (see also Figs. |5(d)| and |5(b)[ ). 




-0.04 



50 100 

Bt/h 



150 



(b) Correlation structure factor S^ 1 ^ vs. rota- 
tional time for 21 sites with 77 = 5, 10. 



L = 21,7] = 10,/3 DC = 1.9 




10" 



0.5 1 1.5 

huj/B 



(d) Squared modulus of Fourier transform of 
Si 10 -' vs. rotationally scaled frequency for L = 21 
sites and 77 = 10. Many new frequencies appear, 
in particular the Rabi frequency and double the 
Rabi frequency, denoted with arrows. 



Figure 5. Dependence of the asymptotic behavior of structure factors on the lattice height 77. 



state (Fig. 4(a)), and the emergent time scale for reaching this equilibrium is shorter than it was for 
77 = 1 by roughly a factor of four. As the lattice height is then increased to rj — 10, the populations 
return to the trend of r] = 1, again converging to quarter filling with a time scale comparable to that of 
77 = 1 (Fig. |4(c)| ). This illustrates the fact that the emergent time scale r is not, in general, a monotonic 
function of the parameters of the lattice. 

While the dynamics of the site-averaged rotational state populations are superficially similar for 
77 = 1 and 77 = 10, the underlying physics is not identical, as can be seen by comparing Figs. 2(f)[ 
4(b)) andgdj These fi gures display the squared modulus of the Fourier transform of the site-averaged 



number in the J = state. The only significant frequency observed for 77 = 1 is the Rabi frequency 
Q ~ 0.064:B/h. In contrast, the 77 = 5 case has numerous other characteristic frequencies. As we raise 
the lattice height to 77 = 10, the frequencies that arose for r] = 5 remain, even though the overall visual 
trend of the site-averaged number reflects that of the single-frequency rj — 1 behavior. While we do not 
explicitly see the new frequencies in the site-averaged number, we do see them in the structure factors. 
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T 

L 


V 


tB/Ti 


Asymp. S.E. 




Asymp. b.b. 


9 


1 


A 1 A f\ A 

414.04 


0.72% 


O r\ o A 

398.4 


n r 1 oy 

0.51%) 


9 


2 


224.32 


1.79/0 


149.9 


1.36%) 


9 


3 


117.5 


1.86% 


126.7 


1.03%) 


9 


10 


613.00 


1.07% 


1079.66 


14.09% 


10 


1 


259.96 


0.76% 


240 


rv a f a 07 

0.6454% 


10 


4 


140.70 


1.19% 


72.04 


0.60% 


10 


10 


526.21 


0.88% 


396.46 


1.018% 


21 


-i 

i 


756.18 


3.13% 


110.68 


0.96/0 


21 


5 


177.53 


1.62% 


75.18 


0.902% 


21 


10 


716.21 


2.96% 


244.09 


2.82% 



Table 3. Emergent time scales r and tq and their fit asymptotic standard errors for various lattice 
heights and system sizes. 



An example is Fig. |5(b) , which clearly displays the 2Vt frequency behavior of the correlation structure 
factor S^ 01 ) for rj = 10. This frequency, which we easily pick out in the site-averaged number's Fourier 
transform, can also be seen in the Fourier transform of see Fig. |5(d) , 



We find that the emergent time scale r does not depend strongly on the size of the system L, even 
though the distribution of molecules on the lattice is, in general, quite different for different numbers 
of sites, as can be seen by comparing Figs. 2(a) and 2(e), Examining Fig. 2(c) and Table 03 the L = 10 



case has a smaller r than either of the odd L cases. We think this has to do with the filling being 
exactly 1/2 and not, strictly speaking, with the number of lattice sites, as the L = 9 and L = 21 cases 



have fillings less than 1/2. We see this clearly by comparing Fig. [4] with Figs. 2(a), 2(c), and 2(e) 



Fig. [4] displays (noO)/AT, a quantity which is independent of filling but dependent, in general, on the 
number of lattice sites. There is a weak dependence on the number of lattice sites. On the other hand, 
Figs. 2(a)[ 2(c), and 2(e) display (hw)/L, a quantity which is independent of the number of lattice sites 
but dependent, in general, on the filling. There is a marked difference between L = 10, which has filling 
of 5/10 = 1/2 and the others, which have fillings< 1/2, but there is not a significant difference between 
L = 9 and L = 21, which have fillings of 4/9 and 10/21, respectively. 

The dependence of r on the filling is also evidenced by the correlation structure factor S^ 01 ) in 
Fig. 3(f)[ which shows that there is a stronger correlation between the J = and J = 1 states for 
exactly half filling than for fillings less than half, regardless of the system size. Half filling is known 
to be important in the extended Bose Hubbard model, where it marks the introduction of the charge 
density wave phase. We thus interpret this greater correlation structure factor as the appearance of a 
dynamic charge density wave phase between rotational states at half filling. 

This is in contrast to the usual behavior, where the structure factors S^°°) and are nonzero 
whenever there is nonzero occupation of the particular rotational state and the structure factor 
is much smaller-essentially zero, see Figs. 3(a)| and 3(e), These results for the structure factors means 
that the J = and J = 1 states tend to lie on top of one another, and not to "checkerboard" 
with a different rotational state occupying alternating sites. This is due to the fact that the Rabi 
flopping time scale is much shorter than the dipole-dipole time scale, meaning that the population cycles 
before there is sufficient time for the molecules to rearrange to a configuration which is energetically 
favorable with respect to the dipole-dipole term. However, because the population in each rotational 
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level asymptotically reaches some nonzero value, we do see a small amount of rearrangement after many 
Rabi periods for any filling, corresponding to a nonzero S^ 01 \ Note that this rearrangement does not 
affect the site-averaged numbers, but rather the distribution of rotational states among the lattice sites. 
This asymptotic distribution emerges on time scales longer than we have considered, and is more prone 
to finite size effects than the site-averaged quantities, so we do not make a conjecture about it here. 
We find that the Q-measure saturates as 

Q = Qmax - AQe~ U/TQ , (65) 

with a different time scale tq, see Fig. [Hand Table [31 We also find that the saturation time scale of the 
Q-measure is not, in general, a monotonic function of the lattice height 77, as shown in Fig. [H 




80 100 

Bt/h 

(a) Dependence of the population damping time scale r 
on the number of lattice sites. When we remove the 
dependence on the filling by dividing through by the 
total number, we see that there is little difference in the 
time scales with which systems of different size approach 
dynamic equilibrium. Contrast Figs. |2(a)[ |2(c) , and 2(e)[ 
which display a profound dependence on filling when the 
dependence on lattice sites has been removed. 




(b) Dependence of spatial entanglement on number of 
lattice sites. We see that systems of different size have 
different spatial entanglement in their static ground state. 
The time scale of the Q-measure saturation, tq, is shorter 
for L = 10 than it is for the odd L cases. This follows the 
general trend of r and tq responding correspondingly to 
changes in the Hamiltonian parameters, and so we associate 
this shorter time scale partially with the filling, not entirely 
with the system size. 



Figure 6. Dependence of emergent time scales on number of lattice sites. 



This time scale is different from the time scale r at which the populations approach an asymptotic 
equilibrium, though both time scales respond similarly to changes in the Hamiltonian parameter, see 
Table [31 For example, if tq gets larger as a parameter is changed then r also gets larger, as illustrated 
in Figs. [Hand [7(b) . The time scale tq displays a stronger dependence on the number of lattice sites L 
than r, as can be seen in Figs. [HandED This is because r describes a quantity that has been averaged 
over sites, while tq does not. 



5. Conclusions 



We have presented and derived a novel lattice Hamiltonian, the Molecular Hubbard Hamiltonian 
(MHH). The MHH is a natural Hamiltonian for connecting theoretical studies of the dynamics of 
quantum phase transitions to near-term experimental setups using ultracold molecular gases. We 
presented a case study of this new Hamiltonian for hard core bosonic molecules at half filling. Starting 
from an initial condition of half filling in the J = 0, M = state, we found that initial large 
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L = 10,/3 D c= 1-9 




20 40 60 80 100 120 140 160 



Btfh 

(a) Dependence of spatial entanglement on lattice height. 
Note that the spatial entanglement and its associated time 
scale are not monotonic functions of the lattice height. 
Note also that the entanglement of the static ground state 
appears to be largely insensitive to the lattice height. 



L= 10,/3 D c= 1.9 




Bt/h 

(b) Dependence of the site-averaged number on the lattice 
height. Note that the emergent time scale r is not a 
monotonic function of the lattice height. Note also that r 
responds in the same way that tq does to changes in the 
lattice height. 



Figure 7. Dependence of emergent time scales on lattice height. 

oscillations in the system self-damp to an asymptotic equilibrium which consists of a lattice height 
and filling-dependent spatially entangled superposition of dressed states. This occurs on an emergent 
time scale r which can not be predicted from the single molecule theory. We showed that r depends 
non-monotonically on lattice height, weakly on lattice size, and strongly on filling (as apparent in 
simulations with odd and even numbers of sites). We also discovered a separate emergent time scale tq 
which describes how quickly the many body spatial entanglement saturates. We demonstrated that tq 
and r respond similarly to changes in the Hamiltonian parameters and that tq depends on the filling, 
the lattice size, and, non-monotonically, on the lattice height. In addition to these emergent time scales, 
we studied the time-dependent structure factors and their frequency-domain Fourier transforms. 

In future studies we will consider different filling factors, DC field strength to rotation ratios 
/?DC? an d initial conditions, as well as polarized and unpolarized spin- 1/2 fermionic molecules. In 
addition, we will use multiscale methods to study how the emergent time scale demonstrated above 
compares to experimental time scales for physical systems, and thereby make quantitative predictions 
for experiments. 

We acknowledge useful discussions with Deborah Jin, Heather Lewandowski, and Jun Ye. This 
work was supported by the National Science Foundation under Grant PHY-0547845 as part of the NSF 
CAREER program. 

Appendix A. Single molecule physics 

Relationship between operators in space-fixed and molecule-fixed coordinate systems 

It is well known that the representation of the angular momentum operators in a molecule-fixed 
coordinate frame lead to the anomalous commutation relations [J i? = —ihcijkJk [48] . The simplest 
way to avoid this trouble is to transform all expressions into the space-fixed frame where the angular 
momentum operators satisfy the normal commutation relations [J^ = ihcijkJk |49j . If the molecule- 
fixed axes are obtained by rotation of the space- fixed axes through the Euler angles {0, 0, x} f28] 
(which we collectively abbreviate as (R)), then the component of a fc t/l -rank spherical tensor T that 
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has projection p along the space- fixed z axis, denoted (T)^\ can be expressed in terms of the molecule 
fixed components as 

(T)L k) =E^W(T)f } , (a.i) 



pq 



where (R)* is the complex conjugate of the pq element of the fc t/l -rank rotation matrix (Wigner 
D-matrix). To avoid confusion, we will label all space- fixed components with the letter p and all 
molecule-fixed components with q. From the orthogonality of the rotation matrices we have the inverse 
relationship 

(T) ( q k) =E^?(R)(r)f (A.2) 

V 

= E(-ir g ^,-,(Rr(rt (a.3) 

V 

Rotational Hamiltonian 

In the rigid rotor approximation the rotational Hamiltonian is simply 

£ rot = B3 2 , (A.4) 
where we have defined the rotational constant B = 1/2/xr^, with ji the molecule's reduced mass and r e 



its equilibrium internuclear separation. Typical values of B are ~ 60/i GHz [50] . This Hamiltonian has 
eigenvalues B J (J + 1) and eigenstates | JM), with J the total angular momentum and M its projection 
along the internuclear axis. 

DC Field Term 

The dipole moment of a polar molecule in a rotational eigenstate is zero in an average sense due to 
the spherical symmetry of the rotational Hamiltonian. We break this symmetry by introducing a DC 
electric field along the space- fixed z axis, with Hamiltonian 

H BC = -d • £ BC , (A.5) 

where £dc is the electric field amplitude. The field defines the spherical space-fixed axis p — 0, and the 
molecule- fixed internuclear axis defines q = 0. We transform between them using a first-rank rotation 
matrix as outlined above: 

£ DC = - (d)^ Sue- (A.6) 

The matrix elements of the DC Hamiltonian in the basis which diagonalize the rotational Hamiltonian 
Eq. ( IA.5D are 

(J',M'\H DC \J,M) = -d£^(2J+l)(2J' + l)(-l) M (A.7) 
x 

where we use the notation (...) for the Wigner 3-j symbol [28]. Note that the symbol d refers to the 
permanent dipole moment of a molecule, and is not to be confused with the dipole operator denoted by 
d. We refer to the basis which simultaneously diagonalizes the Rotational and DC Hamiltonians as the 
"dressed basis," and we denote the kets that span this basis by \£; JM), where the labels J and M are 
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J ) 




1 




—M 





M' ) 


' o 
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the zero field values of the corresponding quantum number and the symbol 8 is a reminder that these 
kets are superpositions of field free rotational states and DC field. 

The effects of the DC field can be clearly seen by considering the dressed state wavefunctions, 
energies, and dipole moments to lowest order in perturbation theory in the dimensionless parameter 
/3 DC = d£ DC /£>, the ratio of the field energy to the rotational level splitting: 



\S, J, M) = | J, M) - _y 17J _ r | J — 1, M) + A 4{J+lf _ 1 \J + 1, Af> , (A.i 



J(J + 1) -3M 



2 



A p(2) _ ^DC 

JM 2B [j(J + l)(2J-l)(2J + 3)J ' 1 J 

,M|d„; JM)jd _ -|£ _ fec ^g^ , (A.10) 

where AEjm is the lowest non-zero shift in the energy. 

The DC field mixes states of different J, breaking the (2 J + l)-fold degeneracy of the rotational 
Hamiltonian, and so J is no longer a good quantum number. In the case of a z-polarized field, M 
remains a good quantum number, and a degeneracy persists for all states with the same \M\. This 
mixing aligns the molecule with the field, inducing a nonzero dipole moment. This means of orienting 
polar molecules, known as "brute force" orientation, works well for molecules that both have a large 
dipole moment and can be efficiently rotationally cooled [SI]. While more effective means of orienting 
molecules using intense laser fields are known [52], they complicate the theoretical discussion and the 
experimental setup, and so we do not consider them here. 

In larger fields the rotational levels become deeply mixed, which allows states that are weak-field 
seeking in low fields to become high-field seeking in high fields [53]. The actual mixing of rotational 
levels vs. /?dc is depicted in Fig. lA2l for the lowest three dressed levels. We note that there always exists 
a field Sr such that the lowest R dressed states' dipole moments are all positive, as this is important 
to ensure the stability of a collection of dipoles. The universal curve of the induced dipole moments 



(in units of d) vs. /?dc of the first two dressed rotational manifolds are shown in Figure 1(a), The 



universal curve of the dressed state energies energies (in units of B) vs. /?dc is shown in Figure l(b)[ 

1 \T 

For reference, /?dc = 1 corresponds to a field of roug hly 1.93^ for B - 60ft GHz and d - 9 D. 

Expanding the field operators in Eq. (pQ) in a Wannier basis of dressed states centered at a particular 
discrete position as described in Eq. (Ill), we find 

j 

H rot + #00 = ^2 ^J^M^S.JM i (A. 11) 

J M=-J 



where E JM is the energy of the \£; J, M) dressed state (see Fig. 2(a)) and tis^jm is the number operator 
associated with this same state. 

If the DC field were aligned at a small angle 9 a to the z field of the trap (say, in the xz plane), 
then small dipole moments mixing M 1 — M ± 1 states would arise and the M' = M dipoles would 
decrease slightly (we can view them as being in an effective field of £ g q = cos9 a £^c)- Treating the new 
contribution perturbatively in the small parameter sin# a /? D c, we find the lowest order couplings to the 
ground state 

. - . s'm9 a d£ / 49sin 2 £L 9 \ /A . 

(£; 00\H DC \£; 1 ± 1) ~ —±- ( 1 - — X cj , (A.12) 
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/5dc 

(a) Scaled induced dipole moments vs. scaled DC field energy. 



4 




/5dc 



(b) Scaled dressed energies vs. scaled DC field energy. 



Figure Al. Dressed state dipole moments and energies. Note that the J = 1,M = resonant 
dipole moment changes from weak- field seeking to high- field seeking at /fee ~ 5. All rotational states 
have a field where this transition occurs, and the dipole tends monotonically towards unity after this 
field. The (10|d|00) dipole moment (and all transition dipole moments, generically) tends towards zero 
monotonically as /fee increases. Note also that the energetic differences between rotational levels are 
smallest at zero field and grow monotonically thereafter. 




^ 0.6 





(£; 10|00)|' 
(8] 10 10) ' 
(£; 10 20) ' 






2 4 6 8 10 

(a) Composition of 1 st dressed (b) Composition of 2 nd dressed (c) Composition of 3 rd dressed 
state. state. state. 

Figure A2. Compositions of dressed states vs. scaled rotational energy. The states become deeply 
mixed in large fields, and that the dressed state \£ ; JM) whose zero field value is | JM) does not always 
have the greatest overlap with | JM) for all /fee • The field strength where the first dressed state changes 
from weak-field to high- field seeking, /fee = 5, is also roughly the place where its overlap with the 1 00) 
field- free level is greater than the overlaps with all other field-free levels. 



and associated timescale TQ a for occupation of M^0 states from the ground state, 

V6h V6 h 



sm9 a d£ 



49 sin" 1 



1440 



/? DC sin 9 a B 



(A.13) 



AC Field Term 

An AC microwave field of frequency cj resonantly drives transitions between two DC dressed 
states \£\ J'M') and |£; JM) with energy difference [Eji M > — E JM ) /h ~ w provided the induced dipole 
moment (8 ; J'M'\d\£] JM) is nonzero. Two states separated by an energy difference AE that is off- 
resonant from the driving field (i.e. AE 3> uj) will also be coupled, albeit much more weakly. In our 
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system we resonantly couple the lowest two dressed rotational levels, |£; 10) and |£;00). We consider 
the case of z polarization, in which the effective Hamiltonian in the dressed Wannier basis is 

H A c (t) = -7r sin (cut) ftjM (4 ; j,m«£; J+i,m + h.c) , ( A. 14) 



JM 



where 



ft JM = S AC (S] J, M|d|£; J + 1, M)/% . (A.15) 

is the Rabi frequency. This is the frequency with which the populations of a two-level system cycles. 
In experiments, the AC field has spatial curvature on the order of cm which is negligible on the /im 
system size scale. 

In the absence of couplings between sites, the physics of the system is determined by the on-site, 
single-molecule physics. The percentage population of each component in both the |£; J, M) dressed 
and \ JM) field-free bases are shown below for one Rabi period. In these plots only the \£; 10) and 
\8] 00) dressed states are considered, which is close to the actual behavior when all other states are far 
off-resonant. Each site undergoes Rabi flopping independently of the others. Figs. |3(a) and |3(b)| show 
this behavior for /?dc — 1.900 and (3ac = d^Ac/B = 0.200, giving a Rabi period of 27r/fi o = S6.5H/B. 




(a) Populations of the dressed states vs. rotational time. The 
small amplitude rapid oscillations occur on the time scale 
l/o;, and are often averaged away via the rotating wave 
approximation. The large amplitude oscillations occuring on 
the time scale l/^oo that periodically transfer the population 
between \£; 00) and \S; 10) are the characteristic "Rabi 
oscillations" of a driven two-level system. 




(b) Populations of the field-free states vs. rotational time. 
The 1 20) state is occupied because both |£;00) and \S; 10) 
have a nonzero projection with this state due to the mixing 
from the DC field, see Fig. IA21 It is apparent from 



comparison with Fig. 3(a) that the dressed basis greatly 
simplifies the AC term in the Hamiltonian. 



Figure A3. Resonant AC field induced population cycling in the dressed and field-free bases. 



Appendix B. Convergence 

Single Molecule Considerations 

Each dressed state \£ ; J, M) is, in principle, an infinite linear combination of field free states 

oo 

\£;J,M)= 52cj,\J',M). (B.l) 

J'=0 

Numerically, we must have a finite upper bound to the sum in Eq. ( IB. II ). which we call J cut . This does 
not cause difficulty in practice, as the overlap of a dressed state \£\ JM) with a field-free state \ J'M) 
diminishes rapidly as J f differs more greatly from J. We find the coefficients in Eq. ( IB. II ). as well as 
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the dressed state energies and dipole moments by simultaneously diagonalizing the rotational and DC 
field Hamiltonians in a basis consisting of the first J cut rotational levels. Because TEBD scales poorly 
with the on-site dimension, we form as small an on-site basis as possible by keeping the eigenvectors 
corresponding to the R lowest dressed levels. To form a proper basis, we must renormalize these 
eigenvectors (which, for z-polarized field, does not change their orthogonality). We now demonstrate 
the convergence of these two procedures 

To show convergence of the first procedure, we plot the difference between the energy of the J th 
rotational state calculated for a particular value of J cut = i and one higher value, AEj (i) as a function of 
i. The results for various field strengths are shown in Figures 1(a)} 1(b) . We see very fast convergence for 
the low fields (e.g. /?dc — 1-9) of interest. In our numerics we use J cut = 25, which ensures convergence 
for any of the /?dc considered. 



/3 D c= 1.9 



{3dc = 20 




10 12 14 16 




8 10 12 14 16 



Figure Bl. Convergence with respect to DC dressing rotational state cutoff. As few as 7 field- free levels 
are needed for the weak field /3d c = 1-9 to have the dressed state energies of interest converge to machine 
precision (left panel), and even a large DC field /?dc = 20 requires only 12 field- free levels for the energy 
to converge (right panel). 



To determine convergence with respect to the second procedure, examine Figs. [2(a)p(b)| which 



show 



p ( j R) = i - j: m jo\io}? 



R-l 



(B.21 



the amount of the total dressed wave function norm \ (S] J0\S] J0)\ 2 that lies outside of the first R field- 
free rotational levels for R = 3 and R = 4, respectively. For R — 4 the renormalization of the first three 
rotational levels is a very small effect for the /?dc we consider, and the fourth level is not populated to 



any appreciable extent during time evolution for any /? D c (see Fig ]3(a) ), so we expect that keeping the 
R — 4 lowest levels will give sufficient accuracy. By direct simulation, we find six digit accuracy in the 
suite of quantum measures defined in Sec. 13.21 specifically, we compare R = 3 to R = 4. 



Many Body Considerations 

There are also convergence issues that are inherent to the TEBD algorithm. The first, called the 
Schmidt error, is the error that arises from truncating the Hilbert space at each time step. We can 
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Figure B2. Convergence with respect to local dimension cutoff. Dressed states with greater J lose 
more of their norm in truncation, as mixing occurs most strongly with adjacent J. Also, as the field is 
increased, the states become more deeply mixed, and so all states lose more of their norm. Truncating 
the local basis at the J = 3 dressed level incurs at most a 1% loss of norm for any of the states that are 
appreciably populated during time evolution (right panel). 



parameterize the error per step in terms of the entanglement cutoff parameter x as 

rf = l- t W) 2 (B-3) 

a t =l 

where is a vector containing the eigenvalues of the reduced density matrix obtained by tracing over 
all sites but /, and ol\ is the local index that entangles the site / with the rest of the system, with smaller 
oil states having greater weight. We find that, among the measures we use, the one that is the most 
sensitive to x is the Q-measure, which we plot for four values of x i n Fig- IB3L Increasing x improves 
the accuracy over longer times, but there is always a time after which the measure begins to deviate. 
This is the normalization drift alluded to in Sec. 13.11 The ^-dependent time after which the Schmidt 
error dominates is referred to as the runaway time [Bj. In the case study of Sec. we used x = 50 
for all simulations, which gives the Q-measure accurately to within four decimal places over the time 
scales considered. 

The second intrinsic source of error in TEBD is due to the Trotter-Suzuki expansion of the 
propagator [37]. We parameterize this error in terms of St, the time step. When we halve the time step 
from that used in the simulations above (= 27t/(133cj)), we find no change in the measures to the ninth 
digit. It is clear that the Schmidt error discussed above is the chief source of error in our simulations. 

To extract the emergent time scales defined in Eqs. f |63l ) and f |65l ). we used two different methods. 
The first is the nonlinear curve fitting routine "fit" in gnuplot. The second is the "Nonlinear Regression" 
package in Mathematica 6.0. Both methods use nonlinear regression, which fits the data to a specified 
nonlinear function of the model parameters. The goodness of the fit is quantified by the asymptotic 
standard errors of the model parameters, which gives the standard deviation of each parameter. A 
low percent asymptotic error means that the model parameters cannot be adjusted very far without 
noticeably changing the goodness-of-fit. Both gnuplot and Mathematica returned the same values for 
the emergent time scales to within the stated asymptotic standard error. 
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Figure B3. Convergence with respect to entanglement cutoff parameter. The left figure shows the 
spatial entanglement measure Q for various values of the TEBD entanglement cutoff parameter \. As 
X is increased, Q remains close to its true value for longer. In the right figure we plot the log of the 
absolute difference in Q for two values of x divided by its arithmetic mean. We see at least four-digit 
accuracy for the largest values of x we consider. Note also that even small values of x are accurate for 
short times. 
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